---
title: "IMCE and VIMP estimation for Yemen"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
#library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc) # write_html()
library(cjbart)
library(Cairo)

dat <- read_csv("data/cj_data_Yemen.csv")
dat <- na.exclude(dat)

dat <- dat |>
  dplyr::select(selected, Activity, Armament,
                Patron, Service, respondentIndex,
                RecordedDate,
                Q1_1, Q1_3, 
                Q2_1, Q2_2, Q2_3, Q2_5, Q2_7,
                Q4_1, Q4_2, Q4_3,
                Gender, D6, D2,
                Q4, Q6,
                Q4_1, Q4_2, Q4_3)
str(dat)
dat[,c(2:5)] <- data.frame(lapply(dat[,c(2:5)],as.factor))

# covariate for VIMP
dat <- dat |> 
  mutate(islamism = Q1_1,
         tribalism = Q1_3,
         patronege = Q4,
         us = Q2_1,
         iran = Q2_2,
         saudi = Q2_3,
         uae = Q2_5,
         axis = Q2_7,
         huthi = Q4_1,
         GPC = Q4_2,
         GPC_IRGC = Q4_3,
         piety = Q6)


# Ceasefire: unexpected event during survey
#dat$RecordedDate <- as.Date(stri_match_first_regex(dat$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat$ceasefire <- NA
dat$ceasefire[dat$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat$ceasefire[as.Date("2025-01-17") <= dat$RecordedDate] <- "Post"
dat$ceasefire <- factor(dat$ceasefire)

# Gender
dat <- dat |> 
  mutate(Gender = case_when(Gender == 1 ~ "Male",
                            Gender == 2 ~ "Female"))
dat$Gender <- as.factor(dat$Gender)

# Age group
dat <- dat |> 
  mutate(Age_group = case_when(D2 == 1 ~ "18-20",
                               D2 == 2 ~ "21-30",
                               D2 == 3 ~ "31-40",
                               D2 == 4 ~ "41-50",
                               D2 == 5 ~ "51-60"))
dat$Age_group <- as.factor(dat$Age_group)

# Sect
dat <- dat |> 
  mutate(Sect = case_when(D6 == 1 ~ "Sunni",
                          D6 == 2 ~ "Zaydi",
                          D6 == 3 ~ "Others"))
dat$Sect <- as.factor(dat$Sect)

# Islamism
dat <- dat |> 
  mutate(Islamism = case_when(Q1_1 >= 67 ~ "High",
                              Q1_1 <= 66 & Q1_1 >= 34 ~ "Middle",
                              Q1_1 <= 33 ~ "Low"))
dat$Islamism <- as.factor(dat$Islamism)


# Tribalism
dat <- dat |> 
  mutate(Tribalism = case_when(Q1_3 >= 67 ~ "High",
                               Q1_3 <= 66 & Q1_3 >= 34 ~ "Middle",
                               Q1_3 <= 33 ~ "Low"))
dat$Tribalism <- as.factor(dat$Tribalism)


# US
dat <- dat |> 
  mutate(US = case_when(Q2_1 >= 4 ~ "Like",
                        Q2_1 == 3 ~ "Neither",
                        Q2_1 <= 2 ~ "Dislike"))
dat$US <- as.factor(dat$US)

# Iran
dat <- dat |> 
  mutate(Iran_subgroup = case_when(Q2_2 >= 4 ~ "Like",
                          Q2_2 == 3 ~ "Neither",
                          Q2_2 <= 2 ~ "Dislike"))
dat$Iran_subgroup <- as.factor(dat$Iran_subgroup)

# Saudi
dat <- dat |> 
  mutate(Saudi = case_when(Q2_3 >= 4 ~ "Like",
                           Q2_3 == 3 ~ "Neither",
                           Q2_3 <= 2 ~ "Dislike"))
dat$Saudi <- as.factor(dat$Saudi)

# UAE
dat <- dat |> 
  mutate(UAE = case_when(Q2_5 >= 4 ~ "Like",
                         Q2_5 == 3 ~ "Neither",
                         Q2_5 <= 2 ~ "Dislike"))
dat$UAE <- as.factor(dat$UAE)


# Axis
dat <- dat |> 
  mutate(Axis = case_when(Q2_7 >= 4 ~ "Like",
                          Q2_7 == 3 ~ "Neither",
                          Q2_7 <= 2 ~ "Dislike"))
dat$Axis <- as.factor(dat$Axis)

# Houthi Allied parties supporter
dat <- dat |> 
  mutate(Houthiallied_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Houthiallied_party_supporter = ifelse(Houthiallied_party_supporter >= 7, "Houthis Allied Party Supporter", "Others"))
dat$Houthiallied_party_supporter <- as.factor(dat$Houthiallied_party_supporter)

# Patronage
dat <- dat |> 
  mutate(Patronage = case_when(Q4 >= 3 ~ "High",
                               Q4 <= 2 ~ "Low"))
dat$Patronage <- as.factor(dat$Patronage)

# High patronage and Houthi allied supporter
dat <- dat |>
  mutate(Patron_Houthi = ifelse(Houthiallied_party_supporter == "Houthis Allied Party Supporter" & Patronage == "High", "Patron_Houthis", "Other"))
dat$Patron_Houthi <- as.factor(dat$Patron_Houthi)


# Piety
dat <- dat |> 
  mutate(Piety = case_when(Q6 <= 4 ~ "Low",
                           Q6 >= 5 & Q6 <= 6 ~ "Middle",
                           Q6 >=7 ~ "High"))
dat$Piety <- as.factor(dat$Piety)

```


# IMCE
```{r, fig.height = 8, dpi = 300, warning = FALSE}
set.seed(89)

# modeling
cj_model <- cjbart(data = dat,
                   Y = "selected",
                   type = "rating",
                   id = "respondentIndex")

# IMCE estimation
het_effect <- IMCE(data = dat,
                   model = cj_model,
                   attribs = c("Activity", "Armament", "Patron","Service"),
                   ref_levels = c("Armed resistance", "Disarmament", "Yemeni government","Distribution of daily goods"),
                   cores = 32,
                   method = "bayes",
                   keep_omce = FALSE)
summary(het_effect)
plot(het_effect)


p1 <- plot(het_effect,
     covar = c("iran"),
     plot_levels = c("Military enforcement", "Iran","Political party")) +
  labs(title = "Yemen")
p1
ggsave(filename = "result/IMCE_Yemen.png", 
       plot = p1, 
       dpi = 600, 
       width = 7, 
       height = 4)

plot(het_effect,
     covar = c("Iran_subgroup"),
     plot_levels = c("Military enforcement", "Iran","Political party"))

plot(het_effect,
     covar = c("ceasefire"),
     plot_levels = c("Military enforcement","Iran","Reconstruction of infrastracture"))
     
```


# VIMP
```{r, fig.height = 7, fig.width = 8, dpi = 300, warning = FALSE}
# random forest model
vimp <- het_vimp(het_effect,
                 covars = c("iran",
                            "Iran_subgroup",
                            "ceasefire"),
                 cores = 32)

plot(vimp,
     covars = c("iran",
                "Iran_subgroup",
                "ceasefire"))

# plot heatmap
covar_lookup <- data.frame(covar_name = c("iran",
                                          "Iran_subgroup",
                                          "ceasefire"),
                           covar_label = c("Iran",
                                          "Iran_subgroup",
                                          "Ceasefire"))

final_results <- vimp$results %>% 
  left_join(covar_lookup, by = c("covar" = "covar_name")) %>% 
  mutate(lab_text = ifelse(importance < 10,"",round(importance,2)))

p20 <- ggplot(final_results, aes(x = covar_label, y = Level, fill = importance)) +
  facet_grid(Attribute ~ ., space = "free", scales = "free", switch = "y") +
  geom_tile() +
  scale_fill_gradient(low="white", high="firebrick2") +
  labs(x = "Subject covariate", y = "Attribute-level", fill = "Importance",
       title = "VIMP of Yemen") +
  theme(text = element_text(size = 14),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
p20

ggsave(filename = "result/VIMP_Yemen.pdf", 
       plot = p20, 
       dpi = 300, 
       width = 5, 
       height = 6)
ggsave(filename = "result/VIMP_Yemen.png", 
       plot = p20, 
       dpi = 600, 
       width = 5, 
       height = 6)
```
